/*---------------------------------------------------------------------------
	ZaRan	-	A Totallly Automatic CFD Software
	Copyright (C) ,Since 2020
-------------------------------------------------------------------------------
License
	This file is part of ZaRan.

!	@file		gird.h
!	@brief	the purpose of this file.
!	@author	Chen Jie.
\*---------------------------------------------------------------------------*/
#include"../include/vanleer.h"
#include"Basic/include/CommonPara.h"
VarList VanLeer::Plus(vector<double*>& localPrimitive, const vector<double>& coord_trans_coef)
{
	double& density = *localPrimitive[FieldDataIndex::density];
	double& u = *localPrimitive[FieldDataIndex::u];
	double& v = *localPrimitive[FieldDataIndex::v];
	double& w = *localPrimitive[FieldDataIndex::w];
	double& p = *localPrimitive[FieldDataIndex::p];
	double c = sqrt(gama * p / density);
	double coef_x = coord_trans_coef[0];
	double coef_y = coord_trans_coef[1];
	double coef_z = coord_trans_coef[2];
	double coef_t = coord_trans_coef[3];
	double d_coef = sqrt(coef_x * coef_x + coef_y * coef_y + coef_z * coef_z);
	double ub = u * coef_x + v * coef_y + w * coef_z + coef_t;
	ub /= d_coef;
	coef_x /= d_coef;
	coef_y /= d_coef;
	coef_z /= d_coef;
	double M = ub / c;
	VarList flux_plus;
	flux_plus.SetVarNum(5);
	if (M <= -1.0)
	{
		flux_plus.SetToSameNumber(0);
	}
	else if (M >= 1.0)
	{
		double E = 0.5 * density * (u * u + v * v + w * w) + p / (gama - 1);
		flux_plus.SetVariable(0, density * ub);
		flux_plus.SetVariable(1, density * ub * u + p * coef_x);
		flux_plus.SetVariable(2, density * ub * v + p * coef_y);
		flux_plus.SetVariable(3, density * ub * w + p * coef_z);
		flux_plus.SetVariable(4, ub * (E + p) - p * coef_t);
	}
	else
	{
		double FM = 0.25 * density * c * (M + 1) * (M + 1);
		flux_plus.SetVariable(0, FM);
		flux_plus.SetVariable(1, FM * (coef_x * (-ub + 2 * c) / gama + u));
		flux_plus.SetVariable(2, FM * (coef_y * (-ub + 2 * c) / gama + v));
		flux_plus.SetVariable(3, FM * (coef_z * (-ub + 2 * c) / gama + w));
		flux_plus.SetVariable(4, FM * (ub * (-ub + 2 * c) / (gama + 1) + 2 * c * c / (gama * gama - 1) + 0.5 * (u * u + v * v + w * w) - coef_t * (-ub + 2 * c) / gama));
	}
	flux_plus *= d_coef;
	return flux_plus;
}


VarList VanLeer::Minus(vector<double*>& localPrimitive, const vector<double>& coord_trans_coef)
{
	double& density = *localPrimitive[FieldDataIndex::density];
	double& u = *localPrimitive[FieldDataIndex::u];
	double& v = *localPrimitive[FieldDataIndex::v];
	double& w = *localPrimitive[FieldDataIndex::w];
	double& p = *localPrimitive[FieldDataIndex::p];
	double c = sqrt(gama * p / density);
	double coef_x = coord_trans_coef[0];
	double coef_y = coord_trans_coef[1];
	double coef_z = coord_trans_coef[2];
	double coef_t = coord_trans_coef[3];
	double d_coef = sqrt(coef_x * coef_x + coef_y * coef_y + coef_z * coef_z);
	double ub = u * coef_x + v * coef_y + w * coef_z + coef_t;
	ub /= d_coef;
	coef_x /= d_coef;
	coef_y /= d_coef;
	coef_z /= d_coef;
	double M = ub / c;
	VarList flux_minus;
	flux_minus.SetVarNum(5);
	if (M >= 1)
	{
		flux_minus.SetToSameNumber(0);
	}
	else if (M <= -1)
	{
		double E = 0.5 * density * (u * u + v * v + w * w) + p / (gama - 1);
		flux_minus.SetVariable(0, density * ub);
		flux_minus.SetVariable(1, density * ub * u + p * coef_x);
		flux_minus.SetVariable(2, density * ub * v + p * coef_y);
		flux_minus.SetVariable(3, density * ub * w + p * coef_z);
		flux_minus.SetVariable(4, ub * (E + p) - p * coef_t);
	}
	else
	{
		double FM = -0.25 * density * c * (M - 1) * (M - 1);
		flux_minus.SetVariable(0, FM);
		flux_minus.SetVariable(1, FM * (coef_x * (-ub - 2 * c) / gama + u));
		flux_minus.SetVariable(2, FM * (coef_y * (-ub - 2 * c) / gama + v));
		flux_minus.SetVariable(3, FM * (coef_z * (-ub - 2 * c) / gama + w));
		flux_minus.SetVariable(4, FM * (ub * (-ub - 2 * c) / (gama + 1) + 2 * c * c / (gama * gama - 1) + 0.5 * (u * u + v * v + w * w) - coef_t * (-ub - 2 * c) / gama));
	}
	flux_minus *= d_coef;
	return flux_minus;
}

VarList VanLeer::UpdateFlux(vector<double*>& localPrimitive, vector<vector<double*>>& neiborPrimitiveVec, CoordTrans& localCoordTrans)
{
	VarList  Flr, Fcl, Fcr, Frl;
	VarList  Glr, Gcl, Gcr, Grl;
	VarList  Hlr, Hcl, Hcr, Hrl;
	vector<double>coord_trans_coef(4);
	coord_trans_coef[0] = localCoordTrans.xi_x();
	coord_trans_coef[1] = localCoordTrans.xi_y();
	coord_trans_coef[2] = localCoordTrans.xi_z();
	coord_trans_coef[3] = localCoordTrans.xi_t();
	Flr = Plus(neiborPrimitiveVec[2], coord_trans_coef);
	Fcl = Minus(localPrimitive, coord_trans_coef);
	Fcr = Plus(localPrimitive, coord_trans_coef);
	Frl = Minus(neiborPrimitiveVec[0], coord_trans_coef);

	coord_trans_coef[0] = localCoordTrans.eta_x();
	coord_trans_coef[1] = localCoordTrans.eta_y();
	coord_trans_coef[2] = localCoordTrans.eta_z();
	coord_trans_coef[3] = localCoordTrans.eta_t();
	Glr = Plus(neiborPrimitiveVec[3], coord_trans_coef);
	Gcl = Minus(localPrimitive, coord_trans_coef);
	Gcr = Plus(localPrimitive, coord_trans_coef);
	Grl = Minus(neiborPrimitiveVec[1], coord_trans_coef);
	if (localCoordTrans.dim() == 2)
	{
		Hlr.SetVarNum(5);
		Hlr.SetToSameNumber(0);
		Hcl = Hcr = Hrl = Hlr;
	}
	else
	{
		coord_trans_coef[0] = localCoordTrans.zeta_x();
		coord_trans_coef[1] = localCoordTrans.zeta_y();
		coord_trans_coef[2] = localCoordTrans.zeta_z();
		coord_trans_coef[3] = localCoordTrans.zeta_t();
		Hlr = Plus(neiborPrimitiveVec[5], coord_trans_coef);
		Hcl = Minus(localPrimitive, coord_trans_coef);
		Hcr = Plus(localPrimitive, coord_trans_coef);
		Hrl = Minus(neiborPrimitiveVec[4], coord_trans_coef);
	}

	VarList F, G, H;
	F = Fcr - Flr + Frl - Fcl;
	G = Gcr - Glr + Grl - Gcl;
	H = Hcr - Hlr + Hrl - Hcl;
	return Fcr - Flr + Frl - Fcl + Gcr - Glr + Grl - Gcl + Hcr - Hlr + Hrl - Hcl;

}
VarList VanLeer::UpdateFlux(vector<double*>& localPrimitive, vector<vector<double*>>& neiborPrimitiveVec)
{

	VarList  Flr, Fcl, Fcr, Frl;
	VarList  Glr, Gcl, Gcr, Grl;
	VarList  Hlr, Hcl, Hcr, Hrl;
	Flr = Split(neiborPrimitiveVec[2], 1, 1.4);
	Fcl = Split(localPrimitive, 2, 1.4);
	Fcr = Split(localPrimitive, 1, 1.4);
	Frl = Split(neiborPrimitiveVec[0], 2, 1.4);

	Glr = Split(neiborPrimitiveVec[3], 1, 1.4);
	Gcl = Split(localPrimitive, 2, 1.4);
	Gcr = Split(localPrimitive, 1, 1.4);
	Grl = Split(neiborPrimitiveVec[1], 2, 1.4);

	Hlr = Split(neiborPrimitiveVec[5], 1, 1.4);
	Hcl = Split(localPrimitive, 2, 1.4);
	Hcr = Split(localPrimitive, 1, 1.4);
	Hrl = Split(neiborPrimitiveVec[4], 2, 1.4);
	return Fcr - Flr + Frl - Fcl + Gcr - Glr + Grl - Gcl + Hcr - Hlr + Hrl - Hcl;
}
VarList VanLeer::Split(vector<double*>& localPrimitive, double key, double gama)
{
	double& density = *localPrimitive[FieldDataIndex::density];
	double& u = *localPrimitive[FieldDataIndex::u];
	double& v = *localPrimitive[FieldDataIndex::v];
	double& w = *localPrimitive[FieldDataIndex::w];
	double& p = *localPrimitive[FieldDataIndex::p];
	double c = sqrt(gama * p / density);
	double xMa = u / c;
	double yMa = v / c;
	double zMa = w / c;
	double E = 0.5 * density * (u * u + v * v + w * w) + p / (gama - 1);
	VarList flux;
	flux.SetVarNum(5);
	if (key == 1)//X+
	{

		if (xMa <= -1.0)
		{
			flux.SetToSameNumber(0);
		}
		else if (xMa >= 1.0)
		{
			flux.SetVariable(0, density * u);
			flux.SetVariable(1, density * u * u + p);
			flux.SetVariable(2, density * u * v);
			flux.SetVariable(3, density * u * w);
			flux.SetVariable(4, u * (E + p) + u * p);
		}
		else
		{
			double FM = 0.25 * density * c * (xMa + 1) * (xMa + 1);
			flux.SetVariable(0, FM);
			flux.SetVariable(1, FM * ((-u + 2 * c) / gama + u));
			flux.SetVariable(2, FM * v);
			flux.SetVariable(3, FM * w);
			flux.SetVariable(4, FM * (((1 - gama) * u * u + 2 * (gama - 1) * u * c + 2 * c * c) / (gama * gama - 1) + 0.5 * (u * u + v * v + w * w) + u * (-u + 2 * c) / gama));
		}
	}
	else if (key == 2)
	{
		if (xMa <= -1.0)
		{
			flux.SetVariable(0, density * u);
			flux.SetVariable(1, density * u * u + p);
			flux.SetVariable(2, density * u * v);
			flux.SetVariable(3, density * u * w);
			flux.SetVariable(4, u * (E + p) + u * p);

		}
		else if (xMa >= 1.0)
		{
			flux.SetToSameNumber(0);
		}
		else
		{
			double FM = -0.25 * density * c * (xMa - 1) * (xMa - 1);
			flux.SetVariable(0, FM);
			flux.SetVariable(1, FM * ((-u - 2 * c) / gama + u));
			flux.SetVariable(2, FM * v);
			flux.SetVariable(3, FM * w);
			flux.SetVariable(4, FM * (((1 - gama) * u * u - 2 * (gama - 1) * u * c + 2 * c * c) / (gama * gama - 1) + 0.5 * (u * u + v * v + w * w) + u * (-u - 2 * c) / gama));
		}
	}
	else if (key == 3)
	{
		if (yMa <= -1.0)
		{
			flux.SetToSameNumber(0);
		}
		else if (yMa >= 1.0)
		{
			flux.SetVariable(0, density * v);
			flux.SetVariable(1, density * v * u);
			flux.SetVariable(2, density * v * v + p);
			flux.SetVariable(3, density * v * w);
			flux.SetVariable(4, v * (E + p) + v * p);
		}
		else
		{
			double FM = 0.25 * density * c * (yMa + 1) * (yMa + 1);
			flux.SetVariable(0, FM);
			flux.SetVariable(1, FM * u);
			flux.SetVariable(2, FM * ((-u + 2 * c) / gama + u));
			flux.SetVariable(3, FM * w);
			flux.SetVariable(4, FM * (((1 - gama) * v * v + 2 * (gama - 1) * v * c + 2 * c * c) / (gama * gama - 1) + 0.5 * (u * u + v * v + w * w) + v * (-v + 2 * c) / gama));
		}
	}
	else if (key == 4)
	{

		if (yMa <= -1.0)
		{
			flux.SetVariable(0, density * v);
			flux.SetVariable(1, density * v * u);
			flux.SetVariable(2, density * v * v + p);
			flux.SetVariable(3, density * v * w);
			flux.SetVariable(4, v * (E + p) + v * p);
		}
		else if (yMa >= 1.0)
		{
			flux.SetToSameNumber(0);
		}
		else
		{
			double FM = -0.25 * density * c * (yMa - 1) * (yMa - 1);
			flux.SetVariable(0, FM);
			flux.SetVariable(1, FM * u);
			flux.SetVariable(2, FM * ((-v - 2 * c) / gama + v));
			flux.SetVariable(3, FM * w);
			flux.SetVariable(4, FM * (((1 - gama) * v * v - 2 * (gama - 1) * v * c + 2 * c * c) / (gama * gama - 1) + 0.5 * (u * u + v * v + w * w) + v * (-v - 2 * c) / gama));
		}
	}
	else if (key == 5)
	{
		if (zMa <= -1.0)
		{
			flux.SetToSameNumber(0);
		}
		else if (zMa >= 1.0)
		{
			flux.SetVariable(0, density * w);
			flux.SetVariable(1, density * w * u);
			flux.SetVariable(2, density * w * v);
			flux.SetVariable(3, density * w * w + p);
			flux.SetVariable(4, w * (E + p) - w * p);
		}
		else
		{
			double FM = 0.25 * density * c * (zMa + 1) * (zMa + 1);
			flux.SetVariable(0, FM);
			flux.SetVariable(1, FM * u);
			flux.SetVariable(2, FM * v);
			flux.SetVariable(3, FM * ((-w + 2 * c) / gama + w));
			flux.SetVariable(4, FM * (((1 - gama) * w * w + 2 * (gama - 1) * w * c + 2 * c * c) / (gama * gama - 1) + 0.5 * (u * u + v * v + w * w) + w * (-w + 2 * c) / gama));
		}
	}
	else if (key == 6)
	{

		if (zMa <= -1.0)
		{
			flux.SetVariable(0, density * w);
			flux.SetVariable(1, density * w * u);
			flux.SetVariable(2, density * w * v);
			flux.SetVariable(3, density * w * w + p);
			flux.SetVariable(4, w * (E + p) - w * p);

		}
		else if (zMa >= 1.0)
		{
			flux.SetToSameNumber(0);
		}
		else
		{
			double FM = -0.25 * density * c * (zMa - 1) * (zMa - 1);
			flux.SetVariable(0, FM);
			flux.SetVariable(1, FM * u);
			flux.SetVariable(2, FM * v);
			flux.SetVariable(3, FM * ((-w - 2 * c) / gama + w));
			flux.SetVariable(4, FM * (((1 - gama) * w * w - 2 * (gama - 1) * w * c + 2 * c * c) / (gama * gama - 1) + 0.5 * (u * u + v * v + w * w) + w * (-w - 2 * c) / gama));
		}
	}
	return flux;
}

